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Jammed packings of granular materials differ from systems normally described by statistical mechanics in that they are athermal. 
In recent years a statistical mechanics of static granular media has emerged where the thermodynamic temperature is replaced by 
a configurational temperature X which describes how the number of mechanically stable configurations depends on the volume. 
Four different methods have been suggested to measure X. Three of them are computed from properties of the Voronoi volume 
distribution, the fourth takes into account the contact number and the global volume fraction. This paper answers two questions 
using experimental binary disc packings: First we test if the four methods to measure compactivity provide identical results 
when applied to the same dataset. We find that only two of the methods agree quantitatively. Secondly, we test if X is indeed an 
intensive variable; this becomes true only for samples larger than roughly 200 particles. This result is shown to be due to recently 
found correlations between the particle volumes [Zhao et al, Europhys. Lett., 2012, 97, 34004]. 



1 Is there a well defined configurational tem- 
perature? 

Temperature is the concept that helps us to understand how 
the exchange of energy stored in the microscopic degrees of 
freedom follows from the accompanying change of entropy 
of the involved systems. If we coarse-grain our view to the 
j macroscopic degrees of freedom of particulate systems, such 
as foams or granular gases, we can still define effective tem- 
peratures that describe their dynamics 1,2 . This approach de- 
fines these systems as dissipative; kinetic energy of the parti- 
cles is irrecoverably lost to microscopic degrees of freedom. 

In the absence of permanent external driving such sys- 
tem will always evolve towards a complete rest. And in the 
presence of boundary forces or gravity this rest state will 
be characterized by permanent contacts between the particles 
which allow for a mechanical equilibrium. Shahinpoor- and 
Kanatani 4 were the first to suggest that such systems might 
still be amenable to a statistical mechanics treatment. Sam 
Edwards and coworkers 5,6 have then developed this idea into a 
full statistical mechanics of static granular matter by using the 
ensemble of all mechanically stable states as a basis. A crucial 
requirement of such an approach is the existence of some type 
of excitation which lets the system explore the phase space of 
all allowed static configurations. This could e.g. be realized 
by tapping, cyclic shear, or flow pulses of the interstitial liq- 
uid. While there are promising results, the feasibility of this 
approach is still under debate 7 . 

A second key concept in Edwards' approach is the replace- 
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ment of the energy phase space by a volume phase space 
where the volume function W(q) takes the role of the the clas- 
sical Hamilton function. The configurations q represent the 
positions and orientations of all grains. One can then define 
an analog to the partition function Z{X) : 



Z(X) 



-W(q)/X 



6(q)dq 



(1) 



where 0(q) limits the integral to mechanically stable configu- 
rations. X is the configurational temperature or compactivity, 
which is defined as X = dV/dS. The configurational entropy 
S corresponds to the logarithm of the number of mechanically 
stable configurations for a given set of boundary condition. 

S is neither known from first principles (except for model 
systems 8,9 ) nor can it be measured directly. Therefore "ther- 
mometers" measuring X have to exploit other relationships; 
four different ways to do so have been suggested. In this pa- 
per we will test all four of them using the same dataset of 
mechanically stable disc packings. 

First X can be determined from the steady state volume 
fluctuations using an analog to the relationship between spe- 
cific heat and energy fluctuations — ~ ; we will refer to this 
compactivity as Xy p • A second method is based on the prob- 
ability ratio of overlapping volume histograms 14,15 , allowing 
us to compute Xoh- A third way 16,17 computes Xr from 
Gamma function fits to the volume distribution. Finally, it 
has been suggested very recently—, that an analysis based on 
"quadrons" 19 instead of Voronoi cells leads to an expression 
for Xq that involves the average particle volume and the con- 
tact number. 

This difference in approach immediately raises the question 
if these four methods provide identical results when applied 
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Fig. 1 Correlations in binary disc packings, a) experimental im- 
age and b) the corresponding Voronoi tessellation, c) The correla- 
tion between two Voronoi volumes as a function of their distance L 
(measured in units of small disc diameters d s ). Dense packings ex- 
hibit anti-correlations for L larger than approx. 2.5 d s . Reproduced 
from 23 . 



to the same experiment. There have been two previous exper- 
iments addressing this question: McNamarra et al 15 found 
that Xvf and Xoh agree for tapped packings of approxi- 
mately 2000 glass spheres. The same result has been reported 
by Puckett and Daniels 20 for compressed disc packings. 

The full description of a granular packing has to take the 
boundary stresses into accountiSi 2 !. The stress dependence of 
the entropy then gives rise to a tensorial temperature named 
angoricity. Angoricity has been computed from numerics 
using the overlapping histogram method 22 and from experi- 
ments using both fluctuations and overlapping histograms 20 . 
As the experiments described below are performed in an open 
cell with gentle driving, the boundary stresses can be assumed 
to be small and constant; we therefore exclude angoricity from 
our further analysis. 

2 Correlations in the Voronoi-volume of disc 
packings 

The three methods to compute Xvf, Xoh, and Xr all start 
from the Voronoi-volume distribution. However, it has been 
shown that the Voronoi volumes inside a sample are corre- 



lated 24 . We recently measured the spatial extension of this 
correlations and demonstrated the existence of additional anti- 
correlations between the volumes 23 . These correlations raise 
the question if X is indeed an intensive parameter i.e. if its 
value is independent of the number N of particles analyzed? 
As the Voronoi- volumes reported in reference 23 will be the 
basis of this study we quickly recap the relevant experimen- 
tal procedures and results, more details can be found in the 
original publication. 



2.1 Experimental setup 

The experiment is performed in an air-fluidized bed filled with 
a binary mixture of Teflon discs with d s = 6mm and d\ = 9 
mm diameter. Data sets of 8000 different mechanically stable 
configurations are prepared by repeated air pulses. Changing 
pressure and duration of the air pulses allows us to control the 
average packing fractions </> in the ranges of 0.818 to 0.838. 
After each flow pulse the discs come to a complete rest and 
are then imaged with a CCD camera (Fig.[T]a). After detect- 
ing the particle centers, the Voronoi volume V of each disc 
is determined (Fig. H]b). We then compute the free volume 
Vf = V — Vmi n for each particle with v m i n being the volume 
the grain would occupy in a hexagonal packing of identical 
discs. This step allows us to superimpose the results for small 
and large discs in the subsequent analysis. 

Our results will also depend on the packing fraction of the 
loosest possible packing </)rlp =0.81 1. This value is averaged 
over 10 packings prepared by slowly settling the discs in a 
manually decreased air flow. 



2.2 Correlations in binary disc packings 

The correlation between the Voronoi volumes can be mea- 
sured using: 



corr(L) 



(( v f,i- v f,i)( v fJ ~ v fd)) 



(2) 



Here i,j are two points in the packing which separated by a 
distance L. The free volumes at these points are Vfj and Vfj\ 
(...) indicates averages over all 8000 different packings and 
240 pairs of points z, j within each packing, a'j is the variance 
of the free volumes averaged over the two points. 

Figure [T]c) depicts corr(L) for two different packing frac- 
tions. At low values of (j) only positive correlation be- 
tween Voronoi cells are found. Above 4>ac = 0.8277 anti- 
correlations appear for L larger than approximately 2.5 small 
particle diameters. We will show below that these (anti-) cor- 
relations control how X becomes intensive. 
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3 Compactivity X VF measured from volume 
fluctuations 

This methods starts from the assumption of a Boltzmann like 
probability distribution: 



-W(cO/X 

p(q,*)= z{x) Q(q) 



(3) 



From equation [3] follows for the probability to observe a 
certain volume V at a given compactivity X: 



p(V,X) = y 5(W(q) - V>(q,X)dq 



= V{V) 



-v/x 



(4) 



Z{X) 



where the density of states V(V) = J 5(W(q) - V)9(q)dq 
counts the number of the mechanical stable states available at 
volume V. Using equation |4] we can determine the average 
volume V(X) as 

p-v/x 



V(X) = V{V) 



Z(X) 



VdV 



(5) 



Taking the derivative of eq. |5]with respect to 1/X shows 
that 

3V 



-(y-vy = -*i 



(6) 



3(1/ X) v ' ' ' v 

on the other hand 3V/3(1/X) can be rewritten as 
— X 2 3V/3X from which follows: 



X 



,dV 



(7) 



dx~ uv 

Nowak et ah were the first to suggest that integrating eq. [7] 
provides a way to compute X— : 

rV RL p $y 
~Z2~ 



1 



X 



VF 



(8) 



v 



where the compactivity at the loosest possible packing has 
been set to infinity: X\fllp = oo. First measurements^^ of 
Xy p showed that it indeed measures a material property as it 
depends on the roughness of the particles. It has to be pointed 
out that this method has the epistemological status of a cal- 
culation rule and not of a test of the Boltzmann assumption 
made in eq. [3] This is different for the overlapping histogram 
method described below. 

As described in section 2 we are interested in the evolution 
of X with the size of the analyzed region, in the following re- 
ferred to as cluster. Therefore we continue our study by using 
the normalized average volume per particle v = ^/N^V/V g 
where the sum goes over all N particles inside the cluster, V 
is the Voronoi volume of the individual particles and V g the 
volume occupied by the particle itself. As a consequence of 
this choice also X is dimensionless. 
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Fig. 2 Compactivity Xvf measured from volume fluctuations, a): 
Average volume variance g% versus average volume v measured for 
clusters of size N — 10, 150. The solid curves are power law fits 
which are then used to numerically integrate eq. |8] (b): The com- 
pactivity Xvf computed from eq.[8]for three different cluster sizes, 
(c): The evolution of Xvf with cluster sizes at different values of 0. 
The radius R of the analyzed cluster is proportional to y/N/<j). 



3.1 Experimental results 

Figure [2(a) shows that for our tapped disc packings the vari- 
ance of the volume fluctuation g\ increases monotonically 
with the average volume v (bar meaning the average over all 
8000 taps). To compute Xvf we perform a power law fit to <r~ 
and use the result to integrate equation [8] numerically. Figure 
[2] (b) shows how the resulting Xy p depends on the packing 
fraction at cluster sizes iV of 1, 30, and 150 discs. It is obvi- 
ous that is not an intensive variable at this range of N. 

3.2 Correlations make Xvf non-intensive in small sys- 
tems 

For a more detailed analysis we have plotted in figure fig. EJc) 
the dependence of Xy f on the cluster radius R measured in 
small disc diameters d s . Three features of Xy f become ap- 
parent: 

1. Xvf is growing monotonously for R < 3d s for all val- 
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ues of 0, 

2. for low to intermediate values of < 
plateau, and 



Xy p then reaches a 



3. for the highest packing fraction Xy F first decreases 
slightly before entering a plateau. 

All three points can be understood by considering the influ- 
ence of the volume correlations shown in figure [T] Equation [8] 
computes Xy f from the average variance o% inside the clus- 
ter. This variance can be decomposed in the following way: 



N N 



k m 

l N 



(9) 



Here o\ = (Svk 2 ) is the fluctuation of a single Voronoi cell 
and Xlm^/c (d v kdv m ) := is the volume correlation between 
disc k and all other discs inside the cluster. If k is in the center 
of the cluster, then Ik is proportional to the area enclosed by 
corr(L), the zero line, and L — R in figure \Hc). 

This consideration allows us to explain feature 1 : while N 
goes from 2 to approximately 25 all Ik are positive and grow- 
ing. Consequentially, the variance of the cluster a? is larger 
than the sum of the variance of the individual Voronoi cells 
a k anc * Xy F grows monotonously with N. 

The second feature stems from the fact, that once disk k 
is more than four to five d s away from the boundary of the 
cluster, Ik becomes independent of the cluster size. As the 
relative importance of "boundary discs" decreases with cluster 
size, o\ goes to a constant and Xy f becomes independent of 
N. 

Finally, the slight decrease in Xyp at high packing frac- 
tions and N values between approximately 30 and 150 can be 
attributed to the anti-correlations that appear for <j) larger than 
0.8277; these will decrease Ik slightly again before it reaches 
its plateau value. 

4 Compactivity Xqh measured from overlap- 
ping histograms 

This way to compute compactivity has been first described by 
Dean and Lefevre 14 . It uses pairs of experiments with slightly 
different values of (j) respectively Xoh and computes then the 
ratio of the probabilities to observe the same local volume V. 
If the assumption of a Boltzmann-like probability distribution, 
as expressed in equation!?! holds, this ratio should be exponen- 
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Fig. 3 Compactivity measured with the overlapping histogram 
method, (a) The probability to observe an average volume per par- 
ticle in clusters with 150 disks. The average packing fraction corre- 
sponds to 0.8336 for the red curve, and 0.8353 for the green curve, 
(b) The logarithm of the probabilities to observe a given volume at 
two different compactivities respectively packing fractions is a linear 
function of the volume, which is in accordance with eq. [TT] Cluster 
size is 150 discs, the dashed lines are linear fits, (c) The % 2 values 
provide a goodness of fit test for both linear (green diamonds) and 
parabolic (red triangle) fits to the probability ratios shown in panel 
b). The average \ 2 is 13% smaller for the linear fit, indicating that a 
Boltzmann-like distribution is a better assumption than a Gaussian. 



tial in V: 



P(V,X! 



P(V,X 2 



T>(V)e- v/x i 

V(V)e- v / x 2 
Z(X 2 ) 



Z(X 2 ) 

z(x 1 ) ( 



(10) 



By taking the logarithm on both sides we obtain: 



Q :=ln 



PxJYl 
Px 2 (V) 



1 

X~2 



1 

X~i 



V + ln 



^x 2 



(11) 



Therefore the difference between two compactivities can be 
computed from a line fit of Q versus V as it was first demon- 
strated in 15 . 

Figure Oa) shows the distribution of average volumes for 
two experiments with a packing fraction difference of 0.0017. 
Figure Ob) demonstrates that the ratio Q is indeed a linear 
function of V, as predicted by eq.Qj] By sweeping the exper- 
imentally accessible range of packing fractions, Xqh can be 
determined from the accumulated compactivity differences up 
to an additive constant X . We determine X by setting Xqh 
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Fig. 4 Comparison of Xoh, Xvf, and Xr computed for single 
discs and clusters of 150 particles. In all three cases the compactivity 
of an individual particle is smaller than that of a larger cluster. Xoh 
and Xvf agree quantitatively, Xr is about an order of magnitude 
smaller. 



for the loosest experimental packing to the value of Xy p at 
this volume fraction. 

The resulting Xoh is shown in figure HI The good quanti- 
tative agreement of Xoh and Xvf is not too surprising given 
that a) both methods are derived from the same probability 
distribution (eq.@]) and b) our determination of X . However, 
the Xoh method provides an additional test of assumption [4] 
as we can compare the quality of a linear fit to Q(V) with 
fit functions derived from alternative probability distribution 
functions 15 . A generic candidate would be a Gaussian distri- 
bution which results in parabolic fit. Figure [2c) demonstrates 
however that a linear fit is superior, adding credibility to a 
Boltzmann-like approach. Also note that canceling the den- 
sity of states in equation [TO] implicitly assumes a weak form 
of ergodicity; if the system explores its phase space different 
at the two values of X we can't eliminate V(V). 



5 Compactivity X r 
distribution fit 



measured from a Gamma 



This third method to compute compactivity has been sug- 
gested by Aste and Di Matteo 17,25 ; it differs from the previous 
approaches that it explicitly determines the density of states 
V(V). Based on the observation 16 that experimental Voronoi 
volume distributions can be well fit by V distributions, the pro- 
pose to replace equation [4] with a rescaled /c-Gamma distribu- 
tion 



P(Yf) = w 



v 4 



(fc-1) 



r(fc) 



exp(— k^J-) 



(12) 
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Fig. 5 Applying the V distribution method, (a) T distribution fits to 
the volume distributions for a single discs and a 150 particle clus- 
ters (0 = 0.8175). (b) Evolution of the normalized shape factor with 
cluster size. 



where Vf is the mean free volume (as defined in section 2), 
and k is the shape factor. They then identify 



Xr 



k 



(13) 



or by making use of the fact that variance cr^ of the T function 
is given by: =V 2 /k they derive 



(14) 



By comparing equation|4]and equation [12] we can also iden- 
tify the density of states: 



goo _ Y - 

z(x) ~ x - 



(fc-i) 



T(k) 



(15) 



Figure [5] (a) demonstrates that our volume distribution is 
well fit by a T distribution. We then determine Xr without 
any additive constant using equation [14] Figure H] show Xr 
in comparison to Xvf and Xoh- While there are qualitative 
similarities as monotonous decrease with (j) and the increase 
with N, both the absolute values and the slope of Xr is quite 
different from the other two methods. 
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The increase of Xr with cluster size points again to the in- 
fluence of the volume correlations. In figure [5] (b) the normal- 
ized shape factor k/N is plotted as a function of N. Similar to 
the evolution of Xyp in figure [2(c), Xr will monotonously 
increase towards a plateau for small values of phi and go 
through a maximum before the plateau for the densest pack- 
ings. 

6 Compactivity Xq measured from quadron 
tessellation 

In a very recent paper Blumenfeld et ah 18 have presented an 
analysis of the statistical mechanics of two-dimensional pack- 
ings based on a quadrons 19 as the building blocks of the tes- 
sellation. An advantage of this choice is that the quadrons 
take by design into account the structural degrees of freedom 
of the individual particles. They then write down the partion 
functions of the volume ensemble Zy , the force ensemble Zp, 
and the total partition function Z, showing in this process that 
Z ^ Z v Zp. Finally, they derive from Z an expression for Xq 
of an isostatic packing: 



(a) 



zN + M 



V 



(16) 



V is the average volume per particle, z the average number 
of contact a particle has, and M is the number of boundary 
forces. 

If Xq is derived from Zy instead of Z one obtains: 



(17) 



In comparing our results with this approach we have to ac- 
knowledge a number of differences. Our experimental pack- 
ings are clearly subject to volume forces due to gravity and 
they are also, as shown in figure [6] (a), hyperstatic. I.e. their 
contact number is larger than what is required to fix all their 
mechanical degrees of freedom. 

Generally, frictional particles only become isostatic at Ran- 
dom Loose Packing and zero pressure 26 . Therefore the only 
direct comparison possible is at our RLP value <\> = 0.81 1, the 
corresponding Xq is shown as black circle in figure |6](c). If 
we assume that we can replace z in equation [T6l with z(<p) and 
use the contact numbers displayed in figure [6] (a), we can also 
compute Xq for a larger range of <p. These results are indi- 
cated as the open squares in figure 0(c). 

As Xq is only computed from the average volume per par- 
ticle, it is insensitive to the correlations described in section 
2. On the other hand there exists a finite size effect if Xq 
is only computed from the volume ensemble, ignoring the 
boundary forces. Figure 0(b) shows how the differences be- 
tween Xq computed from the full and the volume ensemble 
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Fig. 6 Compactivity measured with the quadron tesselation method, 
(a) In all our experiments the contact number per particle is above 
the isostatic value of 3. (b) A comparison of Xq computed from 
the full partiton function (eq. Q2] dashed line) and Xq derived from 
only the volume ensemble (eq. [TTJ circles), (c) A comparison of 
Xvf measured for N = 150 and Xq. The filled circle is computed 
for our Random Loose Packing value which is presumably the only 
isostatic point in our dataset. The open squares assume that equation 
[T6lis also valid for hyperstatic packings. 



vanishes with increasing cluster size and consequentially de- 
creasing contribution of M. 

We acknowledge helpful discussions with Karen Daniels 
and Klaus Kassner. 
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